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Abstract:  In  this  study,  we  present  a  novel  modeling  approaeh  whieh  eombines 
ordinary  differential  equation  (ODE)  modeling  with  logieal  rules  to  simulate  an  arehetype 
bioehemieal  network,  the  human  eoagulation  easeade.  The  model  eonsisted  of  five 
differential  equations  augmented  with  several  logieal  rules  deseribing  regulatory  eonneetions 
between  model  eomponents,  and  unmodeled  interaetions  in  the  network.  This  formulation 
was  more  than  an  order  of  magnitude  smaller  than  eurrent  eoagulation  models,  beeause 
many  of  the  meehanistie  details  of  eoagulation  were  eneoded  as  logieal  rules.  We  estimated 
an  ensemble  of  likely  model  parameters  {N  =  20)  from  in  vitro  extrinsie  eoagulation  data 
sets,  with  and  without  inhibitors,  by  minimizing  the  residual  between  model  simulations 
and  experimental  measurements  using  partiele  swarm  optimization  (PSO).  Eaeh  parameter 
set  in  our  ensemble  eorresponded  to  a  unique  partiele  in  the  PSO.  We  then  validated  the 
model  ensemble  using  thrombin  data  sets  that  were  not  used  during  training.  The  ensemble 
predieted  thrombin  trajeetories  for  eonditions  not  used  for  model  training,  ineluding 
thrombin  generation  for  normal  and  hemophilie  eoagulation  in  the  presenee  of  platelets 
(a  signifieant  unmodeled  eomponent).  We  then  used  flux  analysis  to  understand  how  the 
network  operated  in  a  variety  of  eonditions,  and  global  sensitivity  analysis  to  identify  whieh 
parameters  eontrolled  the  performanee  of  the  network.  Taken  together,  the  hybrid  approaeh 
produeed  a  surprisingly  predietive  model  given  its  small  size,  suggesting  the  proposed 
framework  eould  also  be  used  to  dynamieally  model  other  bioehemieal  networks,  ineluding 
intraeellular  metabolie  networks,  gene  expression  programs  or  potentially  even  eell  free 
metabolie  systems. 
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1.  Introduction 

Developing  mathematieal  models  of  bioehemieal  networks  is  a  signifieant  faeet  of  systems  biology. 
Modeling  approaehes  differ  in  their  degree  of  detail,  where  the  ehoiee  of  approaeh  is  often  determined 
by  prior  knowledge,  or  model  requirements  [1].  Ordinary  differential  equation  (ODE)  models  are 
eommon  tools  for  modeling  bioehemieal  systems  beeause  of  their  ability  to  eapture  dynamies  and  eneode 
meehanism.  However,  ODE  models  typieally  eome  with  diffieult  (or  sometimes  impossible)  parameter 
identifieation  problems.  Eor  example,  Gadkar  et  al.  showed  that  even  with  near-perfeet  information,  it 
was  often  impossible  to  identify  all  the  parameters  in  typieal  signal  transduction  models  [2].  However, 
it  is  not  clear  whether  we  aetually  need  preeise  estimates  for  all  model  parameters.  Bailey  suggested 
more  than  a  deeade  ago,  that  aehieving  qualitative  or  even  quantitative  understanding  of  biologieal 
systems  should  not  require  eomplete  struetural  and  parametrie  knowledge  [3].  Sinee  Bailey’s  eomplex 
biology  with  no  parameters  hypothesis,  Sethna  showed  that  model  performanee  is  typieally  sensitive 
to  only  a  few  parameters,  a  eharaeteristie  seemingly  universal  to  multi-parameter  models  referred  to 
as  sloppiness  [4].  Thus,  reasonable  predietions  may  be  possible,  despite  parametrie  uneertainty,  if  a 
few  eritieal  parameters  are  well-defined.  Eor  example,  Tasseff  et  al,  showed  in  a  model  of  Retinoie 
aeid  (RA)  indueed  differentiation  of  HE-60  eells,  that  eorreet  predietions  were  possible  even  when  75% 
of  the  parameters  were  known  only  to  an  order  of  magnitude  [5].  Perhaps  more  importantly,  ODE 
models  require  signifieant  meehanistie  information,  thereby  limiting  their  utility  in  poorly  understood 
systems,  or  eonversely  explode  in  size  when  eonsidering  multiple  pathways  or  subsystems.  Toward 
this  ehallenge,  logieal  modeling  is  an  emerging  paradigm  that  eneodes  eausal  relationships  between 
model  eomponents  using  quasi-meehanistie  non-linear  transfer  funetions  [6].  Eogieal  models  are  highly 
flexible,  and  despite  their  simplieity,  they  have  eaptured  rieh  behaviors  in  a  variety  of  systems  important 
to  human  health  [7-9].  However,  modeling  eomplex  dynamies  with  logieal  models  is  ehallenging. 
Thus,  there  is  an  unmet  need  for  a  third  approaeh  whieh  eombines  ODEs  and  logieal  models,  where 
ODEs  eould  eneode  meehanistie  information,  while  missing  or  ineomplete  meehanistie  knowledge  ean 
be  approximated  using  a  logieal  approaeh. 

In  this  study,  we  developed  a  hybrid  approaeh  whieh  eombined  ODE  modeling  with  logieal  rules  to 
model  a  well  studied  bioehemieal  network,  the  human  eoagulation  system.  Coagulation  is  an  arehetype 
proteolytie  easeade  involving  both  positive  and  negative  feedbaek  [10-12].  Coagulation  is  mediated  by 
a  family  proteases  in  the  eireulation,  called  faetors  and  a  key  group  of  blood  eells,  called  platelets.  The 
eentral  proeess  in  eoagulation  is  the  eonversion  of  prothrombin  (fll),  an  inactive  eoagulation  faetor,  to 
the  master  protease  thrombin  (Ella).  Thrombin  generation  involves  three  phases,  initiation,  amplifieation 
and  termination  [13,14].  Initiation  requires  a  trigger  event,  for  example  vessel  injury,  whieh  leads  to  the 
aetivation  of  faetor  VII  (EVIIa).  Two  eonverging  pathways,  the  extrinsie  and  intrinsie  easeades,  then 
proeess  and  amplify  this  initial  coagulation  signal.  The  extrinsie  easeade  is  generally  believed  to  be 
the  main  meehanism  of  thrombinogenesis  in  the  blood  [15-17].  Initially,  thrombin  is  produeed  upon 
cleavage  of  prothrombin  by  fluid  phase  activated  faetor  X  (EXa),  whieh  itself  has  been  aetivated  by  Tissue 
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Factor/activated  factor  VII  (TF/FVIIa)  [10].  Picomolar  amounts  of  thrombin  then  activate  the  cofactors 
factors  V  and  VIII  (fV  and  fVIII)  and  platelets,  leading  to  the  formation  of  the  tenase  and  prothrombinase 
complexes  on  activated  platelets.  These  complexes  amplify  the  early  coagulation  signal  by  further 
activating  FXa,  and  directly  converting  prothrombin  to  thrombin.  There  are  several  control  points  in  the 
cascade  that  inhibit  thrombin  formation,  and  eventually  terminate  thrombin  generation.  Tissue  Factor 
Pathway  Inhibitor  (TFPI)  inhibits  FXa  formation  catalyzed  by  TF/FVIIa,  while  antithrombin  III  (ATIII) 
neutralizes  several  of  the  proteases  generated  during  coagulation,  including  thrombin.  Thrombin  itself 
also  inadvertently  plays  a  role  in  its  own  inhibition;  thrombin,  through  interaction  with  thrombomodulin, 
protein  C  and  endothelial  cell  protein  C  receptor  (EPCR),  converts  protein  C  to  activated  protein  C 
(APC)  which  attenuates  the  coagulation  response  by  proteolytic  cleavage  of  fV/FVa  and  fVlll/FVllla. 
Termination  occurs  after  either  prothrombin  is  consumed,  or  thrombin  formation  is  neutralized  by 
inhibitors  such  as  APC  or  ATIII. 

Previous  coagulation  models  have  typically  been  formulated  as  systems  of  nonlinear  ordinary 
differential  equations,  using  mass  action  or  more  complex  kinetics,  to  describe  the  rates  of  biochemical 
conversions  [18-22].  Mechanistic  ODE  coagulation  models  from  our  laboratory  [23,24]  were  built  upon 
the  earlier  studies  of  Jones  and  Mann  [25],  Hockin  et  fl/.  [26],  and  later  Butenas  et  a/.  [27]  who  developed 
and  then  subsequently  refined  highly  mechanistic  coagulation  models.  Other  laboratories  have  also 
expanded  upon  Hockin  et  al.  for  example  by  exploring  the  intrinsic  pathway,  the  role  of  stochastic 
fluctuations  in  coagulation  [28],  and  the  dynamics  of  thrombin  mediated  clot  formation  [29]  and 
fibrinolysis  [30].  Other  aspects  of  coagulation  have  also  been  modeled,  such  as  platelet  biochemistry  [31], 
multi-scale  models  of  clot  formation  [32,33],  and  transport  inside  clots  [34].  However,  these  previous 
studies  were  largely  based  upon  extensive  mechanistic  knowledge.  This  is  possible  because  blood, 
while  enormously  complex,  can  be  systematically  interrogated.  Other  systems,  such  as  intracellular 
signaling  networks,  are  much  more  difficult  to  experimentally  interrogate.  Towards  this  unmet 
need,  we  formulated  a  hybrid  modeling  approach  which  combines  ODEs  and  logical  rules  to  model 
biochemical  processes  for  which  a  complete  mechanistic  understanding  is  missing.  We  tested  this 
approach  by  modeling  the  human  coagulation  cascade.  Others  have  also  constructed  reduced  order 
human  coagulation  models.  Recently,  Papadopoulos  and  co-workers  constructed  a  phenomenological 
mathematical  model  for  thrombin  generation  [35].  Using  four  ordinary  differential  equations  and  six 
parameters,  they  derived  an  expression  describing  the  temporal  evolution  of  thrombin  generation  in  a 
variety  of  cases.  The  reduced  order  Papadopoulos  model  showed  good  agreement  with  experimental 
data,  and  underscored  that  model  reduction  is  possible  even  for  complex  positive  feedback  systems 
like  coagulation.  However,  the  Papadopoulos  model  was  focused  on  thrombin  generation,  and  had  a 
lesser  emphasis  on  the  influence  of  physiological  inhibitors  such  as  ATIII  or  the  protein  C  pathway.  In 
this  study,  we  focused  on  building  a  reduced  order  coagulation  model  that  included  the  physiological 
inhibitors  using  a  hybrid  strategy.  The  hybrid  model  consisted  of  only  five  differential  equations 
augmented  with  several  logical  rules.  Thus,  the  model  was  more  than  an  order  of  magnitude  smaller 
than  comparable  purely  ODE  models  in  the  literature.  We  estimated  the  model  parameters  from  in  vitro 
extrinsic  coagulation  data  sets,  in  the  presence  of  ATIII,  with  and  without  the  protein  C  pathway.  We  then 
compared  the  model  predictions  with  thrombin  data  sets,  for  both  normal  and  hemophilic  coagulation, 
that  were  not  used  for  model  training.  Once  validated,  we  performed  flux  and  sensitivity  analysis  on 
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the  model  to  estimate  whieh  parameters  were  eritical  to  model  performance  in  several  conditions.  The 
reduced  order  hybrid  approach  produced  a  surprisingly  predictive  coagulation  model,  suggesting  this 
framework  could  potentially  be  used  to  model  other  biochemical  networks  important  to  human  health. 


Figure  1.  Schematic  of  the  connectivity  of  the  reduced  order  coagulation  model.  A  trigger 
compound,  e.g.,  TF/FVIIa  initiates  thrombin  production  (Flla)  from  prothrombin  (fll).  Once 
activated,  thrombin  catalyzes  its  own  activation  (amplification  step),  as  well  as  its  own 
inhibition  via  the  conversion  of  protein  C  to  activated  protein  C  (APC).  APC  and  tissue 
factor  pathway  inhibitor  (TFPI)  inhibit  initiation  and  amplification,  while  antithrombin  III 
(ATIII)  directly  inhibits  thrombin.  All  inhibition  steps  and  trigger-induced  initiation  were 
modeled  using  a  rule-based  approach.  Likewise,  the  dependence  of  amplification  on  other 
coagulation  factors  was  also  modeled  using  a  rule-based  approach.  The  abundance  of  the 
highlighted  species  (in  the  dashed  boxes)  was  governed  by  an  ordinary  differential  equation. 
All  other  species  were  assumed  to  be  constant. 


2.  Results 

2.1.  Formulation  of  Reduced  Order  Coagulation  Models 

We  developed  a  reduced  order  extrinsic  coagulation  model  to  test  our  hybrid  modeling  approach 
(Figure  1).  The  core  of  our  model  was  based  upon  the  earlier  work  of  Ismagilov  and  coworkers  [36-39], 
where  we  added  initiation,  factor  dependence,  and  specific  inhibition  terms  to  the  earlier  simplified 
model.  A  trigger  event  initiates  thrombin  formation  (Flla)  from  prothrombin  (fll)  through  a  lumped 
initiation  step.  This  step  loosely  represents  the  initial  activation  of  thrombin  by  activated  FXa.  Once 
activated,  thrombin  catalyzes  its  own  formation  (amplification  step),  and  inhibition  via  the  conversion  of 
protein  C  to  activated  protein  C  (APC).  Antithrombin  III  (ATIII)  inhibits  amplification,  while  APC  and 
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tissue  factor  pathway  inhibitor  (TFPI)  potentially  inhibit  both  initiation  and  amplification.  All  initiation 
and  inhibition  processes,  as  well  as  the  dependence  of  amplification  upon  other  coagulation  factors, 
was  approximated  using  our  rule-based  approach  (Figure  2).  Individual  regulatory  contributions  to  the 
activity  of  pathway  enzymes  were  integrated  into  control  coefficients  (u’s)  using  an  integration  rule 
(min/max).  These  control  coefficients  then  modified  the  rates  of  model  processes  at  each  time  step. 
Hill-like  transfer  functions  0  <  f  (Z)  <  I  quantified  the  contribution  of  components  upon  a  target 
process.  Components  were  either  individual  inhibitor  or  activator  levels  or  some  function  of  levels,  e.g., 
the  product  of  factor  levels.  In  this  study,  Z  corresponded  to  the  abundance  of  individual  inhibitors 
or  activators,  with  the  exception  of  the  dependence  of  amplification  upon  specific  coagulation  factors 
(modeled  as  the  product  of  factors).  When  a  process  was  potentially  sensitive  to  multiple  inputs,  logical 
integration  rules  were  used  to  select  which  transfer  functions  influenced  the  process  at  any  given  time. 
In  our  proof  of  concept  model,  we  used  a  winner  takes  all  strategy;  the  maximum  or  minimum  transfer 
function  was  selected  at  any  given  time  step.  However,  other  integration  rules  are  certainly  possible. 
Taken  together,  while  the  reduced  order  coagulation  model  encodes  significant  biological  complexity,  it 
is  highly  compact  (consisting  of  only  five  differential  equations).  Thus,  it  will  serve  as  an  excellent  proof 
of  principle  example  to  study  the  reduction  of  a  highly  complex  human  subsystem. 

2.2.  Identification  of  Model  Parameters  Using  Particle  Swarm  Optimization 

A  critical  challenge  for  any  dynamic  model  is  the  estimation  of  kinetic  parameters.  We  estimated 
kinetic  and  control  parameters  simultaneously  from  eight  in  vitro  time-series  coagulation  data  sets 
with  and  without  the  protein  C  pathway.  The  residual  between  model  simulations  and  experimental 
measurements  was  minimized  using  particle  swarm  optimization  (PSO).  A  population  of  particles 
{N  =  20)  was  initialized  with  randomized  kinetic  and  control  parameters  and  allowed  to  search  for 
parameter  vectors  that  minimized  the  residual.  However,  not  all  parameters  were  varied  simultaneously. 
We  partitioned  the  parameter  estimation  problem  into  two  subproblems  based  upon  the  biological 
organization  of  the  training  data;  (i)  estimation  of  parameters  associated  with  thrombin  formation  in  the 
absence  of  the  protein  C  pathway  and  (ii)  estimation  of  parameters  associated  with  the  protein  C  pathway. 
Only  those  parameters  associated  with  each  subproblem  were  varied  during  the  optimization  procedure 
for  that  subproblem,  e.g.,  thrombin  parameters  were  not  varied  during  the  protein  C  subproblem.  The 
PSO  procedure  was  run  for  20  generations  for  each  subproblem,  where  each  generation  was  1200 
iterations.  The  best  particle  from  each  generation  was  used  to  generate  the  particle  population  for  the 
next  generation.  We  rotated  the  subproblems,  starting  with  subproblem  1  in  the  first  generation. 

The  experimental  training  data  for  parameter  estimation  was  reproduced  from  the  experiments  of 
Butenas  and  co-workers  [40].  In  these  experiments  thrombin  generation  was  initiated  by  FVIIa-TF 
using  mean  plasma  concentrations  of  coagulation  proteins  and  inhibitors.  To  prepare  FVIIa-TF,  TF 
(0.5  nmol/L)  was  relipidated  into  400  /rmol/L  of  phospholipid  vesicles  (PCPS)  by  incubation  in 
20  mmol/L  HEPES,  150  mmol/E  NaCl,  and  2  mmol/E  CaCl2  pH  7.4  (HBS/Ca^+)  for  30  min  at  37  °C. 
The  relipidated  TE  was  incubated  with  10  pmol/E  factor  Vila  for  20  min  to  allow  the  formation  of 
EVIIa-TE.  Eactors  V,  VIII  and  thrombomodulin  (Tm)  (when  protein  C  activation  is  required)  were  added 
to  EVIIa-TE  complex.  Thrombin  generation  was  then  initiated  by  adding  equal  volumes  of  this  mixture 
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with  a  mixture  containing  prothrombin,  factor  IX  and  factor  X,  TFPI,  AT- III  and  protein  C  (added  when 
required),  protein  S  (added  when  required)  and  factor  XI  (added  when  required).  In  the  training  data, 
5  pmol/L  FVIIa-TF  was  used  along  with  200  /rmol/L  of  phospholipid  vesicles  (PCPS)  to  initiate 
thrombin  generation.  All  other  the  coagulation  factors  and  inhibitors  i.e.,  factors  X,  IX,  V,  and 
VIII,  prothrombin,  TFPI  and  AT-III,  protein  C  and  protein  S  (when  applicable)  were  at  their  mean 
plasma  concentrations. 
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Figure  2.  Schematic  of  rule -based  effective  control  laws.  Traditional  enzyme  kinetic 
expressions,  e.g.,  Michaelis-Menten  or  multiple  saturation  kinetics  are  multiplied  by  an 
enzyme  activity  control  variable  0  <  Uj  <  1.  Control  variables  are  functions  of 
many  possible  regulatory  factors  encoded  by  arbitrary  transfer  functions  of  the  form 
0  <  fj  [Z)  <  1.  At  each  simulation  time  step,  the  Vj  variables  are  calculated  by  evaluating 
integration  rules  such  as  the  max  or  min  of  the  set  of  transfer  functions  /i, . . . ,  /„  influencing 
the  activity  of  enzyme  Ej. 
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Figure  3.  Reduced  order  coagulation  model  training  simulations.  Reduced  order  coagulation 
model  parameters  were  estimated  using  particle  swarm  optimization  (PSO)  without  the 
protein  C  pathway  as  a  function  of  prothrombin.  Solid  lines  denote  the  simulated  mean  value 
of  the  thrombin  profile  for  N  =  20  independent  particles,  points  denote  experimental  data. 

The  shaded  region  denotes  the  99%  confidence  estimate  of  the  mean  simulated  thrombin 
value  (uncertainty  in  the  model  simulation).  (A-C)  depict  training  data  and  results  for  150%, 
100%  and  50%  of  physiological  prothrombin  levels  in  the  absence  of  protein  C  pathway.  The 
experimental  training  data  was  reproduced  from  the  study  of  Butenas  et  al.  [40].  Thrombin 
generation  in  these  experiments  was  initiated  using  5  pmol/L  FVIIa-TF  in  the  presence  of 
200  /imol/L  of  phospholipid  vesicles  (PCPS).  As  depicted  in  (A-C)  the  prothrombin  levels 
were  at  150%,  100%  and  50%  of  their  physiological  concentration  in  the  absence  of  protein 
C  pathway.  All  factors  and  control  proteins  in  these  experiments  were  at  their  physiological 
concentration  unless  otherwise  denoted. 
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Figure  4.  Reduced  order  coagulation  model  training  simulations.  Reduced  order  coagulation 
model  parameters  were  estimated  using  particle  swarm  optimization  (PSO)  with  the 
protein  C  pathway  as  a  function  of  prothrombin.  Only  APC  pathway  parameters  were 
allowed  to  vary  in  these  simulations  keeping  the  parameters  estimated  without  protein  C 
pathways  constant.  Solid  lines  denote  the  simulated  mean  value  of  the  thrombin  profile  for 
N  =  20  independent  particles,  points  denote  experimental  data.  The  shaded  region  denotes 
the  99%  confidence  estimate  of  the  mean  simulated  thrombin  value  (uncertainty  in  the 
model  simulation).  (A-C)  depict  training  data  and  results  for  150%,  100%  and  50%  of 
physiological  prothrombin  levels  in  the  presence  of  the  protein  C  pathway.  The  experimental 
training  data  was  reproduced  from  the  study  of  Butenas  et  al.  [40,41].  Thrombin  generation 
in  these  experiments  was  initiated  using  5  pmol/L  FVIIa-TF  in  the  presence  of  200  /rmol/L 
of  phospholipid  vesicles  (PCPS).  As  depicted  in  (A-C)  the  prothrombin  levels  were  at  150%, 
100%  and  50%  of  their  physiological  concentration  in  the  presence  of  protein  C  pathway.  All 
factors  and  control  proteins  in  these  experiments  were  at  their  physiological  concentration 
unless  otherwise  denoted. 


The  reduced  order  coagulation  model  captured  the  role  of  initial  prothrombin  abundance,  and  the 
decay  of  the  thrombin  signal  following  from  ATIII  activity  (Figure  3).  However,  we  systematically 
under-predicted  the  thrombin  peak  and  the  strength  of  ATIII  inhibition  in  this  training  data  set.  On  the 
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other  hand,  with  fixed  thrombin  parameters,  we  eaptured  peak  thrombin  values  and  the  deeay  of  the 
thrombin  signal  (at  least  for  the  150%  HI  ease)  in  the  presenee  of  both  ATIII  and  the  protein  C  pathway 
(Figure  4).  Lastly,  we  were  unable  to  eapture  global  differenees  in  initiation  time  across  separate 
data  sets  with  a  single  ensemble  of  model  parameters.  These  differenees  likely  resulted  from  normal 
experimental  variability.  For  example,  different  thrombin  generation  experiments  within  the  training 
data  (at  the  same  physiologieal  faetor  levels)  had  signifieantly  different  initiation  times  (data  not  shown). 
However,  the  inability  to  globally  eapture  initiation  time  also  highlighted  a  potential  shorteoming  of 
the  initiation  module  within  the  model.  To  eapture  the  variability  in  initiation  time  across  training 
data  sets,  we  ineluded  a  eonstant  time-delay  parameter  (T^)  for  eaeh  data  group.  The  delay  parameter 
was  eonstant  within  a  data  set,  but  allowed  to  vary  across  training  data  sets.  Introduetion  of  the  delay 
parameter  allowed  the  model  to  simulate  multiple  training  data  sets  using  a  single  ensemble  of  model 
parameters.  Taken  together,  the  model  identifieation  results  suggested  that  our  hybrid  approaeh  eould 
reproduee  a  panel  of  thrombin  generation  data  sets  in  the  neighborhood  of  physiologieal  faetor  and 
inhibitor  eoneentrations.  However,  it  was  unelear  whether  the  redueed  order  model  eould  prediet  new 
data,  without  updating  the  model  parameters. 

2.3.  Validation  of  the  Reduced  Order  Coagulation  Model 

We  tested  the  predietive  power  of  the  redueed  order  eoagulation  model  with  validation  data  sets 
not  used  during  model  training.  Two  validation  data  sets  were  used,  thrombin  generation  for  various 
prothrombin  and  ATIII  eoneentrations  with  the  protein  C  pathway,  and  thrombin  generation  in  normal 
versus  hemophilie  plasma  in  the  presenee  of  the  protein  C  pathway.  Lastly,  we  eompared  the  qualitative 
output  of  the  model  to  rFVIIa  addition  in  the  presenee  of  hemophilia.  The  hemophilia  ease  was  an 
espeeially  diffieult  test  as  it  was  taken  from  a  different  study  whieh  used  a  plasma-based  in  vitro  assay 
involving  platelets  instead  of  phospholipid  vesieles  (PCPS).  All  kinetie  and  eontrol  parameters  were 
fixed  for  the  validation  simulations.  The  only  globally  adjustable  parameter  To,  was  fixed  within 
eaeh  validation  data  set  but  allowed  to  vary  between  data  sets.  The  redueed  order  model  predieted 
the  thrombin  generation  profile  for  ratios  of  prothrombin  and  ATIII  in  the  absenee  of  the  protein  C 
pathway  (Figure  5).  Simulations  near  the  physiologieal  range  (fII,ATIII)  =  (100%,  100%)  or  (125%, 
75%)  traeked  the  measured  thrombin  values  (Figure  5B,C).  On  the  other  hand,  predietions  for  faetor 
levels  outside  of  the  physiologieal  range  (fII,ATIII)  =  (50%,  150%)  or  (150%,  50%),  while  qualitatively 
eonsistent  with  measured  thrombin  values,  did  show  signifieant  deviation  from  the  measurements 
(Figure  5A,D).  Likewise,  simulations  of  thrombin  generation  in  normal  versus  hemophilia  (missing  both 
fVIII  and  fIX)  were  eonsistent  with  measured  thrombin  values  (Figure  6).  We  modeled  the  dependenee 
of  thrombin  amplifieation  on  faetor  levels  using  a  produet  rule  {Z  =  fV  x  fX  x  fVIII  x  fIX),  whieh 
was  then  was  integrated  using  a  min  integration  rule  into  the  eontrol  variable  governing  amplifieation. 
Thus,  in  the  absenee  of  fVIII  or  fIX,  the  amplifieation  eontrol  variable  evaluated  to  zero,  and  the  only 
thrombin  produeed  was  from  initiation  (Figure  6B).  However,  the  deeay  of  the  thrombin  signal  was 
underpredieted  in  the  normal  ease  (Figure  6A),  while  the  aetivated  thrombin  level  was  overpredieted 
in  hemophilia  simulations,  although  thrombin  generation  was  far  less  than  normal  (Figure  6B).  Taken 


Thrombin  [nM]  Thrombin  [nM]  Thrombin  [nM] 


Processes  2015,  3 


187 


together,  the  redueed  order  model  performed  well  in  the  physiologieal  range  of  faetors,  even  with 
unmodeled  eomponents  sueh  as  platelet  aetivation  in  the  hemophilia  data  set. 
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Figure  5.  Redueed  order  eoagulation  model  predietions  versus  experimental  data  for  normal 
eoagulation.  The  redueed  order  eoagulation  model  parameter  estimates  were  tested  against 
data  not  used  during  model  training.  Simulations  of  different  levels  of  prothrombin  and  ATIII 
were  eompared  with  experimental  data  in  the  absenee  of  the  protein  C  pathway.  Solid  lines 
denote  the  simulated  mean  value  of  the  thrombin  profile  for  N  =  20  independent  partieles, 
points  denote  experimental  data.  The  shaded  region  denotes  the  99%  eonfidenee  estimate  of 
the  mean  simulated  thrombin  value  (uneertainty  in  the  model  simulation).  (A-D)  predietion 
results  for  (FII,ATIII):  (50%,  150%),  (100%,  100%),  (125%,  75%)  and  (150%,50%)  of 
physiologieal  prothrombin  and  ATIII  levels  in  the  absenee  of  the  protein  C  pathway.  The 
experimental  validation  data  was  reprodueed  from  the  study  of  Butenas  et  al.  [40] .  Thrombin 
generation  in  these  experiments  was  initiated  using  5  pmol/L  FVIIa-TF  in  the  presenee  of 
200  yumol/L  of  phospholipid  vesieles  (PCPS).  As  depieted  in  (A-D)  the  prothrombin  and 
ATIII  levels  were  at  (50%,  150%),  (100%,  100%),  (125%,  75%)  and  (150%,  50%)  of  their 
physiologieal  eoneentrations  in  the  absenee  of  protein  C  pathway.  All  faetors  and  eontrol 
proteins  were  at  their  physiologieal  eoneentration  unless  otherswise  denoted. 
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Figure  6.  Reduced  order  coagulation  model  predictions  versus  experimental  data  with  and 
without  coagulation  factors  VIII  (fVIII)  and  IX  (FIX).  The  reduced  order  coagulation  model 
parameter  estimates  were  tested  against  data  not  used  during  model  training.  Simulations 
of  normal  thrombin  formation  with  ATIII  and  the  protein  C  pathway  were  compared  with 
thrombin  formation  in  the  absence  of  fVIII  and  flX.  Solid  lines  denote  the  simulated  mean 
value  of  the  thrombin  profile  for  N  =  20  independent  particles,  points  denote  experimental 
data.  The  shaded  region  denotes  the  99%  confidence  estimate  of  the  mean  simulated 
thrombin  value  (uncertainty  in  the  model  simulation).  (A,B)  prediction  results  for  normal 
thrombin  generation  and  thrombin  generation  in  hemophilia.  All  factors  and  control  proteins 
were  at  their  physiological  concentration  unless  others  noted.  Coagulation  was  initiated 
with  0.2  nmol/L  FVIIa.  The  experimental  validation  data  was  reproduced  from  the  study  of 
Allen  et  al.  [42]. 
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Figure  7.  Reduced  order  coagulation  model  predictions  of  rFVIIa  administration. 
(A)  Simulations  of  thrombin  formation  in  the  presence  of  ATIII  and  the  protein  C  pathway 
were  conducted  for  a  range  of  trigger  values  (lx-200x  nominal)  in  the  absence  of  fVIII  and 
flX;  (B)  Comparison  of  thrombin  generation  for  normal  versus  hemophilia  for  lOx  nominal 
trigger.  Solid  lines  denote  the  simulated  mean  value  of  the  thrombin  profile  for  N  =  20 
independent  particles.  The  peak  thrombin  time  for  normal  coagulation  (t*)  is  less  than  rFVIIa 
induced  coagulation  in  hemophilia  (t**),  while  the  peak  thrombin  value  was  greater  in  normal 
coagulation.  The  shaded  region  denotes  the  99%  confidence  estimate  of  the  mean  thrombin 
value  (uncertainty  in  the  model  simulation).  All  factors  and  control  proteins  were  at  their 
physiological  concentration  unless  others  noted. 


The  model  ensemble  predicted  a  direct  correlation  between  thrombin  generation  and  rFVIIa 
addition  in  hemophilia  (Figure  7).  In  the  current  model,  we  cannot  distinguish  between  different 
initiation  sources,  e.g.,  TF/FVIIa  versus  rFVIIa,  as  we  have  only  a  single  lumped  initiation 
source  (trigger).  Thus,  we  simulated  the  addition  of  rFVIIa  in  hemophilia  by  removing  fVIII 
and  HX  from  the  model,  and  modulating  the  initial  level  of  trigger.  Simulations  with  a 
baseline  level  of  trigger  were  consistent  with  the  previous  hemophilia  simulations,  where  the 
only  thrombin  produced  was  from  initiation  (Figure  7 A,  lx  trigger).  However,  as  we  increased 
the  trigger  strength,  the  thrombin  profile  began  to  approximate  normal  coagulation,  showing  a 
pronounced  peak  albeit  with  a  slower  peak  time  (Figure  7B,  t**  >  t*).  Further  increases  in 
trigger  strength  resulted  in  decreased  thrombin  peak  time  and  increased  maximum  thrombin  values 
(Figure  7A,  50x  trigger).  Thus,  for  large  trigger  values  (200x  trigger),  the  hemophilic  thrombin  profile 
approximated  normal  coagulation,  where  peak  thrombin  was  achieved  shortly  after  administration  and 
95%  of  the  thrombin  was  gone  by  20  min  after  initiation.  We  performed  flux  analysis  to  understand 
how  the  reduced  order  coagulation  model  balanced  initiation,  amplification  and  termination  of  thrombin 
generation  for  normal  and  hemophilic  coagulation.  Analysis  of  the  reaction  flux  through  the  reduced 
order  network  for  thrombin  generation  in  normal,  hemophilia  and  rFVIIa-treated  hemophilia  identified 
three  distinct  operational  modes  (Figure  8).  We  calculated  the  flux  through  four  lumped  reactions. 
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initiation,  amplification,  thrombin-induced  APC  generation  and  total  thrombin  inhibition  (ineluding  both 
APC  and  ATIII  aetion).  Directly  after  the  addition  of  a  trigger  (e.g.,  TF/FVIIa  or  rFVIIa),  the  lumped 
initiation  flux  was  the  largest  for  all  three  eases.  However,  within  a  few  minutes  enough  thrombin  was 
generated  by  the  initiation  meehanism  to  induee  the  amplifieation  stage.  During  amplifieation,  thrombin 
eatalyzes  its  own  formation  and  inhibition  by  generating  aetivated  protein  C  (APC),  a  potent  inhibitor 
of  the  eoagulation  easeade.  For  normal  eoagulation,  amplifieation  and  thrombin  inhibition  were  the 
dominate  reaetions  by  6  min  after  initiation  (Figure  8,  left).  After  10  min,  the  dominate  reaetion  had 
shifted  to  thrombin  inhibition  (both  ATIII  and  APC  aetion).  In  hemophilia  (missing  both  fVIII  and  fIX), 
the  amplifieation  reaetion  did  not  oeeur,  and  thrombin  was  produeed  only  by  initiation  (Figure  8,  eenter). 
Initiation  was  quiekly  inhibited  by  APC,  and  the  thrombin  level  stabilized  (eventually  deeaying  at  longer 
times  beeause  of  ATIII  activity).  Lastly,  when  50 x  trigger  was  used  to  induee  thrombin  formation  in 
hemophilia  (absenee  of  fVIIFfIX),  initiation  meehanisms  dominated  for  up  to  6  min  following  initiation 
(Figure  8,  right).  Similar  to  hemophilia  alone,  no  amplification  occurred  in  the  50 x  trigger-i-hemophilia 
ease,  and  the  rate  of  thrombin  generation  was  extinguished  by  the  eombined  aetion  of  ATIII  and 
APC.  Taken  together,  the  hybrid  modeling  approaeh  eaptured  the  transition  between  the  modes  of 
thrombin  generation,  as  well  as  the  role  that  inhibitors  play  in  attenuating  the  thrombin  generation 
rate.  Thus,  the  transfer  funetion  approaeh  encoded  the  inhibitory  logie  of  this  easeade  in  the  absenee 
of  speeifie  meehanism. 
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Figure  8.  Reaetion  flux  distribution  as  a  funetion  of  time  for  thrombin  generation  under 
normal  (left),  hemophilia  (eenter)  and  rFVIIa  treated  hemophilia  (right).  Reaetion  flux  was 
ealeulated  for  eaeh  partiele  at  T  =  0, 4,  6,  8, 10, 12, 14  min  after  the  initiation  of  eoagulation. 
Reaetion  fluxes  were  ealeulated  for  eaeh  partiele  in  the  parameter  ensemble  (N  =  20).  Blue 
eolors  denote  low  flux  values  while  red  eolors  denote  high  flux  values. 
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Figure  9.  Global  sensitivity  analysis  of  the  redueed  order  eoagulation  model  with  respeet 
to  the  model  parameters.  (A)  Sensitivity  analysis  of  the  thrombin  peak  time  for  different 
prothrombin  levels  (150%,  100%  and  50%  of  the  physiologieal  value)  as  a  funetion 
of  activated  protein  C;  (B)  Sensitivity  analysis  of  the  thrombin  exposure  for  different 
prothrombin  levels  (150%,  100%  and  50%  of  the  physiological  value)  as  a  function  of 
activated  protein  C.  Points  denote  the  mean  total  sensitivity  value,  while  the  area  around 
each  point  denotes  the  uncertainty  in  the  sensitivity  value.  The  gray  dashed  line  denotes  the 
45°  degree  diagonal,  if  sensitivity  values  are  equal  for  different  conditions  they  will  lie  on  the 
diagonal.  Sensitivity  values  significantly  above  or  below  the  diagonal  indicate  differentially 
important  model  parameters.  The  radius  of  the  shaded  region  around  each  total  sensitivity 
value  was  the  maximum  uncertainty  in  that  value  estimated  by  the  Sobol  method. 


2.4.  Global  Sensitivity  Analysis  of  the  Reduced  Order  Coagulation  Model 

We  conducted  a  global  sensitivity  analysis  to  estimate  which  parameters  controlled  the  performance  of 
the  reduced  order  model.  We  calculated  the  sensitivity  of  the  time  to  maximum  thrombin  (peak  time)  and 
the  thrombin  exposure  (area  under  the  thrombin  curve)  for  different  levels  of  prothrombin,  and  protein 
C  (Figure  9).  Globally,  41%  of  the  parameters  shifted  in  importance  between  the  (fII,PC)  =  (50%,  0%) 
and  (150%,  100%)  cases  for  the  peak  thrombin  time  (Figure  9A).  The  majority  of  these  shifts  involved 
the  interaction  between  increased  prothrombin  and  the  protein  C  pathway,  while  only  5%  were  directly 


Processes  2015,  3 


192 


associated  with  increased  prothrombin  alone.  The  rate  constant  for  thrombin  amplification  was  the 
most  important  parameter  controlling  the  peak  thrombin  time.  While  this  parameter  was  differentially 
important  for  different  prothrombin  levels,  and  in  the  presence  or  absence  of  the  activated  protein  C 
pathway,  it  was  consistently  the  most  sensitive  parameter  in  the  model.  The  saturation  constant  governing 
thrombin  amplification  was  the  second  most  important  parameter,  followed  by  the  initiation  control  gain 
parameter.  Other  important  parameters  influencing  the  thrombin  peak  time  included  the  control  gain  for 
activated  protein  C  formation,  and  the  rate  constant  controlling  ATIII  inhibition  of  thrombin  activity. 
On  the  other  hand,  only  27%  of  the  model  parameters  were  differentially  sensitive  between  the 
(fII,PC)  =  (50%,  0%)  and  (150%,  100%)  cases  for  thrombin  exposure  (Figure  9B).  Of  these  parameters, 
all  of  the  shifts  were  associated  with  the  interplay  between  thrombin  formation  and  the  protein  C 
pathway.  The  rate  constant  controlling  ATIII  inhibition  was  the  most  important  parameter  controlling 
the  thrombin  exposure.  While  this  parameter  was  less  important  in  the  presence  of  protein  C  for 
150%  prothrombin  levels,  it  was  significantly  above  all  other  parameters.  Similar  to  the  peak  time, 
for  150%  prothrombin,  the  control  gain  for  activated  protein  C  formation  was  differentially  important 
along  with  the  rate  constant  controlling  amplification.  However,  the  amplification  parameter  was  much 
less  important  for  thrombin  exposure  vs.  peak  time. 

3.  Discussion 

In  this  study,  we  developed  a  reduced  order  model  of  the  human  coagulation  cascade.  We  modeled 
coagulation  because  it  is  well  studied,  has  a  complex  architecture,  and  has  an  abundance  of  experimental 
data  available  for  model  identification  and  validation.  However,  coagulation  was  just  a  proof  of  concept 
test  of  our  approach.  The  proposed  hybrid  framework  could  also  be  used  to  dynamically  model 
other  biochemical  networks,  including  intracellular  metabolic  networks,  gene  expression  programs 
or  potentially  even  cell  free  metabolic  systems.  The  model  consisted  of  five  differential  equations 
augmented  with  several  logical  rules  describing  regulatory  connections  between  model  components 
and  unmodeled  interactions  in  the  network.  We  estimated  model  parameters  from  in  vitro  extrinsic 
coagulation  data  sets,  in  the  presence  of  ATIII,  with  and  without  the  protein  C  pathway.  To  estimate 
parameters,  the  residual  between  model  simulations  and  experimental  measurements  was  minimized 
using  particle  swarm  optimization  (PSO).  However,  not  all  of  the  model  parameters  were  uniquely 
identifiable,  given  the  training  data.  Instead,  we  estimated  an  ensemble  of  likely  parameter  sets  (N  =  20) 
from  eight  in  vitro  time-series  coagulation  data  sets  with  and  without  the  protein  C  pathway.  Ensemble 
approaches  have  been  used  previously  for  other  signal  transduction  models  [43-47],  and  for  metabolic 
models  [48]  to  estimate  the  impact  of  poorly  constrained  parameter  values  or  poorly  understood  network 
structure  on  simulation  performance.  Thus,  ensemble  approaches  are  common  in  the  dynamic  modeling 
community.  However,  a  unique  feature  of  the  current  study  is  the  direct  connection  between  our  particle 
swarm  approach,  and  the  parameter  ensemble;  each  particle  in  our  swarm  uniquely  corresponded  to 
a  parameter  set  in  our  ensemble.  Thus,  by  constraining  particles  to  operate  in  different  parameter 
regions,  giving  each  particle  a  different  parameter  combination  to  explore,  or  perhaps  even  suppling 
a  different  model  formulation  to  each  particle  we  can  effectively  traverse  through  complex  parameter 
and  model  spaces.  We  validated  the  ensemble  using  thrombin  data  sets  taken  from  multiple  laboratories 
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for  a  variety  of  experimental  eonditions  not  used  during  training.  The  ensemble  predieted  thrombin 
trajeetories  for  eonditions  not  used  for  model  training,  ineluding  thrombin  generation  for  normal  and 
hemophilie  eoagulation  in  the  presenee  of  platelets  (a  signifieant  unmodeled  eomponent).  We  then  used 
flux  analysis  to  understand  how  the  network  operated  in  a  variety  of  eonditions,  and  global  sensitivity 
analysis  to  identify  whieh  parameters  eontrolled  the  performanee  of  the  network.  Flux  analysis  showed 
the  logieal  rules  formulation  eneoded  the  transitions  between  initiation,  amplifieation  and  termination 
of  thrombin  generation.  Sensitivity  analysis  suggested  that  the  amplifieation  rate  eonstant  was  more 
important  to  the  time  to  peak  thrombin,  while  the  ATIII  inhibition  eonstant  eontrolled  thrombin  exposure. 
Taken  together,  the  proposed  hybrid  framework  produeed  a  surprisingly  predietive  model,  suggesting  this 
approaeh  eould  be  used  to  effeetively  model  other  bioehemieal  networks  important  to  human  health. 

Malfunetions  in  eoagulation  ean  have  potentially  fatal  eonsequenees.  Aggressive  elotting  involved 
with  Coronary  Artery  Diseases  (CADs),  eolleetively  aeeounts  for  38%  of  all  deaths  in  North 
Ameriea  [49].  Coagulation  management  during  surgery  ean  also  be  ehallenging,  partieularly  with  the 
inerease  in  elinieal  use  of  antithrombotie  drugs  [50].  Insuffieient  eoagulation  due  to  genetie  disorders 
sueh  as  hemophilia  ean  also  result  in  reeurrent  bleeding.  The  eoagulation  faetors  VIII  (fVIII)  and  IX 
(fIX)  are  defieient  in  Hemophilia  A  and  B,  respeetively  [51-53].  People  with  mild  hemophilia  have 
5%^0%  of  the  normal  elotting  faetor  levels  while  severe  hemophiliaes  have  <1%  [53].  Hemophilia 
ean  be  eontrolled  with  regular  infusions  of  the  defieient  elotting  faetors.  However,  elotting  faetor 
replaeement  sometimes  leads  to  the  formation  of  fVIII  and  fIX  inhibitors  in  vivo  [54].  Alternatively, 
reeombinant  faetor  Vila  (rFVIIa)  has  been  used  to  treat  bleeding  disorders  [55,56]  ineluding  hemophilia 
with  and  without  faetor  VIII/IX  inhibitors  [57].  However,  rFVIIa  requires  frequent  administration 
(every  2-3  h),  and  many  questions  remain  about  its  meehanism  of  aetion,  its  effeetive  dosage  [54], 
and  its  overall  utility  for  the  treatment  trauma-assoeiated  hemorrhage  [58].  In  this  study,  we  did  not 
model  rFVIIa-indueed  eoagulation  direetly.  Rather,  we  modeled  a  general  trigger  whieh  initiated  the 
extrinsie  eoagulation  easeade.  Sinee  we  identified  the  model  using  TF/FVIIa,  inherent  to  our  rFVIIa 
simulations  (and  the  rate  eonstant  governing  initiation)  was  the  presenee  of  TF.  However,  even  with  this 
eomplieation,  the  model  generated  potentially  useful  insight  into  the  rFVIIa  meehanism  of  aetion,  and 
its  possible  shorteomings  espeeially  for  the  treatment  of  hemophilia.  The  addition  of  rFVIIa  direetly 
aetivated  thrombin  through  the  initiation  pathway.  However,  no  amplifieation  of  the  thrombin  signal 
oeeurred  without  fVIII  or  fIX.  Thus,  the  peak  thrombin  signal  was  lower  than  normal  eoagulation,  the 
peak  thrombin  time  was  longer,  and  thrombin  generation  was  eventually  inhibited  by  the  eombined  aetion 
of  ATIII  and  the  protein  C  pathway.  However,  as  the  dose  of  rFVIIa  inereased,  the  peak  thrombin  time 
deereased  (eventually  saturating  around  200 x  nominal  trigger),  and  the  peak  thrombin  value  inereased 
sueh  that  the  thrombin  profile  resembled  normal  eoagulation.  Butenas  et  al.  performed  an  extensive 
in  vitro  study  of  rFVIIa-indueed  thrombin  generation  under  normal  and  hemophilie  eonditions  [59]. 
They  found  qualitatively  similar  trends,  namely  rFVIIa  restored  normal  eoagulation  (even  in  the  absenee 
of  TF)  for  large  enough  rFVIIa  doses,  although  rFVIIa-indueed  eoagulation  in  hemophilia  (even  for 
large  rFVIIa  doses)  lagged  the  normal  profile.  These  results  suggest  that  rFVIIa  administration  alone 
might  not  be  able  to  initiate  normal  eoagulation  in  reeurrent  bleeding,  unless  the  dosage  is  well  above  a 
eritieal  threshold.  However,  defining  this  threshold,  whieh  is  likely  patient  speeifie,  is  diffieult  as  there 
is  tremendous  patient  to  patient  variability  even  with  a  normal  eoagulation  phenotype  [60] . 
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The  hybrid  model  simulations  of  thrombin  formation  were  qualitatively  similar  to  other  published 
coagulation  models.  Many  mathematical  models  of  coagulation  dynamics  have  been  built  upon  the 
Hockin-Mann  model  [26] .  For  example,  Brummel-Ziedins  and  co-workers  incorporated  the  Protein  C 
(PC)  pathway  into  the  original  Hockin-Mann  network  to  investigate  thrombin  generation  in  cases  of 
familial  PC  deficiency  [61].  Simulations  using  this  model  showed  that  PC  mutations  caused  elevated 
thrombin  levels  without  changing  the  initiation  time  or  the  slope  of  thrombin  generation.  This  trend 
was  qualitatively  captured  by  our  reduced  order  model.  For  example,  we  showed  decreased  peak 
thrombin  concentration  in  the  presence  of  PC  pathway  and  similar  thrombin  generation  slopes,  although 
the  initiation  time  was  different  as  these  were  different  experimental  trials.  Danforth  et  al.  simulated 
normal  thrombin  generation  using  5  pM  tissue  factor  with  all  other  factors  at  their  mean  physiological 
level  [60].  The  initiation  time  in  this  simulation  was  approximately  4.4  min.  When  predicting  the 
normal  thrombin  generation  curve  using  0.2  nM  of  TF/FVIIa,  we  showed  an  initiation  time  between 
4-5  min  in  a  platelet  based  system,  and  an  initiation  time  of  3-4  min  in  a  PCPS  system  with  5  pM 
TF/FVIIa,  although  these  times  were  largely  dictated  by  the  time  delay  parameter  TD.  Thus,  while 
were  not  able  to  explicitly  predict  initiation  time,  we  did  bracket  the  initiation  time  predicted  by  the 
Danforth  study.  Mitrophanov  et  al.  showed,  in  a  study  exploring  the  mode  of  action  of  rfVIIa  [62], 
that  increasing  amounts  of  rFVIIa  accelerated  the  maximal  slope  of  thrombin  generation,  and  both  the 
peak  thrombin  and  initiation  times.  While  we  failed  to  capture  the  effect  of  rFVIIa  on  initiation  time,  we 
correctly  predicted  that  changes  in  the  rFVIIa  concentration  affected  the  maximal  thrombin  slope  and  the 
propagation  phase.  Lastly,  detailed  mechanistic  coagulation  models,  for  example  the  model  by  Luan  and 
co-workers  [24]  or  the  recent  model  by  Mitrophanov  et  al.  [30],  often  contain  hundreds  of  proteins  and 
interactions.  Thus,  it  is  unlikely  that  the  reduced  order  hybrid  model  will  replicate  all  the  rich  behavior 
of  these  other  models.  However,  for  qualitatively  different  cases  such  as  normal  versus  hemophilia,  the 
hybrid  approach  gave  similar  results.  For  example,  akin  to  the  hybrid  model,  Luan  et  al.  modeled  the 
normal  and  hemophilia  data  from  Allen  [42].  However,  unlike  the  previous  Luan  et  al.  study,  we  used 
this  data  for  validation  rather  than  training.  Hybrid  model  simulations  of  the  Allen  et  al.  data  set  were 
surprisingly  consistent  with  the  Luan  et  al.  model.  For  example,  the  amplification  of  a  normal  thrombin 
signal  were  comparable,  and  in  hemophilia  both  correctly  predicted  decreased  thrombin  amplification. 
Taken  together,  the  hybrid  model  performance  was  similar  to  other  full  scale  mechanistic  models  in  the 
literature,  although  we  consistently  failed  to  predict  initiation  times  across  data  sets. 

The  performance  of  the  proof  of  principle  coagulation  model  was  impressive  given  its  limited  size. 
However,  there  are  several  issues  that  could  be  further  explored.  First,  the  prediction  of  initiation 
time  should  be  investigated.  We  were  able  to  estimate  initiation  time  within  a  data  set,  but  unable  to 
predict  initiation  time  across  independent  data  sets.  This  suggested  that  we  should  update  the  initiation 
module  to  distinguish  between  different  triggers,  e.g.,  TF/FVIIa  versus  rFVIIa  alone,  and  to  include  key 
biological  milestones  such  as  FXa  activation  (a  prerequisite  to  thrombin  formation).  Next,  there  are 
several  additional  biological  modules  that  could  be  added  to  the  core  model  presented  here.  First,  we 
could  include  thrombin-induced  platelet  activation  and  the  role  of  activated  platelets  in  amplification.  We 
captured  thrombin  generation  data  in  the  presence  of  platelets,  however,  the  initial  shape  of  the  activation 
curve  and  the  time-scale  of  activation  was  not  always  consistent  with  the  data.  Platelets  are  activated  by 
thrombin  through  the  cleavage  of  the  extracellular  domain  of  protease-activated  receptors  (PARs)  on  the 
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platelet  surfaee.  Onee  aetivated,  platelets  play  an  important  role  in  amplifieation,  and  are  key  mediators 
of  the  positive  feedbaek  driving  amplifieation.  Thus,  this  biology  is  a  potentially  important  eomponent 
of  an  expanded  model.  We  should  also  add  the  intrinsie  pathway  to  the  model.  The  intrinsie  pathway  is 
triggered  by  eontaet  aetivation  of  the  plasma  protease  faetor  XI  (fXI)  by  negatively  eharged  surfaees  and 
by  thrombin  and  upstream  faetors  sueh  as  aetivated  plasma  protease  faetor  XII  (FXIIa)  [63,64] .  Aetivated 
platelets  may  also  release  polyphosphate  whieh  direetly  aetivates  fXII  [65].  Arguably  a  minor  player  in 
aeute  bleeding,  eontaet  aetivation  eould  also  be  important  in  other  wound  healing  eontexts.  Finally,  to 
make  the  model  more  elinieally  relevant,  we  should  inelude  the  bioehemieal  proeesses  responsible  for 
elot  formation  and  elot  dissolution  (fibrinolysis).  Clot  formation  is  driven  by  thrombin  aetivity,  while 
fibrinolysis  is  driven  by  plasmin  aetivity,  a  key  enzyme  that  eleaves  fibrin  (one  of  the  main  materials 
in  a  elot).  Similar  to  eoagulation,  fibrinolysis  is  managed  by  several  aetivating  and  inhibitory  faetors 
whieh  eontrol  the  balanee  between  elot  formation  and  dissolution.  Tissue  plasminogen  aetivator  (t-PA) 
and  urokinase  aetivate  plasmin,  along  with  eontaet  pathway  faetors  sueh  as  fXIa.  On  the  other  hand, 
thrombin  aetivatable  fibrinolysis  inhibitor  (TAFI)  inhibits  the  degradation  of  fibrin  by  plasmin.  Also, 
similar  to  eoagulation,  there  is  eonsiderable  fibrinolysis  and  eontaet  pathway  data  sets  that  ean  be  used 
to  train  the  model.  Lastly,  the  ehoiee  of  max/min  integration  rules  or  the  partieular  form  of  the  transfer 
funetions  eould  be  generalized  to  inelude  other  rule  types  and  funetions.  Theoretieally,  an  integration 
rule  is  a  funetion  whose  domain  is  a  set  of  transfer  funetion  inputs,  and  whose  range  is  u  G  [0, 1]. 
Thus,  integration  rules  other  than  max/min  eould  be  used,  sueh  as  the  mean  or  the  produet,  assuming 
the  range  of  the  transfer  funetions  is  always  /  G  [0, 1].  Alternative  integration  rules  sueh  as  the  mean 
might  have  different  properties  whieh  eould  influenee  model  identifieation  or  performanee.  For  example, 
a  mean  integration  rule  would  be  differentiable,  whieh  allows  derivative-based  optimization  approaehes 
to  be  used.  The  partieular  form  of  the  transfer  funetion  eould  also  be  explored.  We  ehoose  a  Hill-like 
funetion  beeause  of  its  prominenee  in  the  systems  and  synthetie  biology  eommunity.  However,  the  only 
mathematieal  requirement  for  a  transfer  funetion  is  that  it  map  a  non-negative  eontinuous  or  eategorieal 
variable  into  the  range  /  G  [0, 1].  Thus,  many  types  of  transfer  funetions  are  possible. 

4.  Materials  and  Methods 

4.1.  Formulation  and  Solution  of  the  Model  Equations 

We  used  ordinary  differential  equations  (ODEs)  to  model  the  time  evolution  of  proteins  {Xi)  in  our 
redueed  order  eoagulation  model: 


(1) 


where  TZ  denotes  the  number  of  reaetions,  fA  denotes  the  number  of  protein  speeies  in  the  model. 


The  quantity  Vj  (x,  e,  k)  denotes  the  rate  of  reaetion  j.  Typieally,  reaetion  j  is  a  non-linear  funetion  of 
bioehemieal  speeies  abundanee,  as  well  as  unknown  kinetie  parameters  k  (/C  x  1).  The  quantity  aij 
denotes  the  stoiehiometrie  eoeffieient  for  speeies  i  in  reaetion  j.  If  aij  >  0,  speeies  i  is  produeed  by 
reaetion  j.  Conversely,  if  (Jij  <  0,  speeies  i  is  eonsumed  by  reaetion  j,  while  =  0  indieates  speeies 
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i  is  not  connected  with  reaction  j.  The  system  material  balances  were  subject  to  the  initial  conditions 
X  (to)  =  Xo,  which  were  specified  by  the  experimental  setup. 

Each  reaction  rate  was  written  as  the  product  of  two  terms,  a  kinetic  term  (fj)  and  a  control  term  (vj) 
that  could  depend  upon  many  regulatory  transfer  functions: 

Tj  (x,  e,  k)  =  fjVj  (2) 

We  used  multiple  saturation  kinetics  to  model  the  reaction  term  ff 


where  denotes  the  maximum  rate  for  reaction  j,  denotes  the  enzyme  abundance  which  catalyzes 
reaction  j,  and  Kjg  denotes  the  saturation  constant  for  species  s  in  reaction  j.  The  product  in 
Equation  (3)  was  carried  out  over  the  set  of  reactants  for  reaction  j  (denoted  as  mj).  The  control 
term  Vj  depended  upon  the  combination  of  factors  which  influenced  the  activity  of  enzyme  i.  Eor  each 
enzyme,  we  used  a  rule-based  approach  to  select  from  competing  control  factors  (Eigure  2).  If  an  enzyme 
was  activated  by  m  metabolites,  we  modeled  this  activation  as: 

Vj  =  max  (fij  (Z) fmj  (Z))  (4) 

where  0  <  fij  (Z)  <  1  was  a  regulatory  transfer  function  that  calculated  the  influence  of  metabolite  i 
on  the  activity  of  enzyme  j.  Conversely,  if  enzyme  activity  was  inhibited  by  m  metabolites,  we  modeled 
this  inhibition  as: 

Vj  =  l-  max  (fij  (Z),...,  f^j  (Z))  (5) 

Eastly,  if  an  enzyme  had  both  m  activating  and  n  inhibitory  factors,  we  modeled  the  control  term  as: 

Vj  =  min  (uj,dj)  (6) 
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Kjs 


(3) 


where: 


Ui 


m^  (/ij  (Z),..  ■  )  fmj  (2)) 

1  -  max  (fij  (Z),...,  fnj  (Z)) 


(7) 

(8) 


The  quantities  and  j~  denoted  the  sets  of  activating  and  inhibitory  factors  for  enzyme  j.  If 
a  process  has  no  modifying  factors,  we  set  Vj  =  1.  There  are  many  possible  functional  forms  for 
0  <  fij  (Z)  <  1.  However,  in  this  study,  each  individual  regulatory  transfer  function  took  the  form: 


fi  (Zj,  kij)  = 


i  +  klz] 


(9) 


where  Zj  denotes  the  abundance  of  the  j  factor  (e.g.,  metabolite  abundance),  and  kij  and  rj  are  control 
parameters,  kij  was  the  species  gain  parameter,  while  rj  was  a  cooperativity  parameter  (similar  to  a 
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Hill  coefficient).  Applying  the  general  framework  to  the  reduced  eoagulation  network  resulted  in  five 
ordinary  differential  equations: 


j,  (j'init'^init  ^amp^amp) 


dxi 

dt 

dX2 

dt 

dX3 

dt 

dxi 

dt 

dx5 

dt 


f'amp'^amp  “f  '^inh,ATIIl'^inh,ATIII 


—  ^  apc^apc 


—  ^  apc'^apc 


—  —T  inh,ATIIlVinh,ATIII 


(10) 

(11) 

(12) 

(13) 

(14) 


where  x  =  (///,  Flla,  PC,  APC,  ATIIIy  .  The  terms  in  the  balanee  equations  denote  eorreeted 
kinetie  expressions  for  initiation,  amplifieation  and  inhibition  proeesses.  The  rate  of  initiation  finu  was 
modeled  as: 


Tinit  =  kinit  {trigger) 


Xi 


k^init,fll  T 


(15) 


where  kinu,  Kinujii  are  the  rate  and  saturation  constants  governing  initiation,  respectively.  The  rate  of 
initiation  was  modified  by  Vinu,  the  eontrol  parameter  governing  initiation.  Initiation  was  sensitive  to  the 
level  of  trigger  (aetivator)  and  TFPI  (inhibitor): 


Vinit  =  min  {TFPI) ,  /+.,  {trigger)) 


(16) 


where  the  transfer  funetions  /  took  the  form  of  Equation  (9).  The  rate  of  thrombin  amplifieation  was 
given  by: 


^amp  kajnp  (3:2) 


Xi 


(17) 


where  fcamp,  Kampjii  denote  the  rate  and  saturation  eonstants  governing  amplifieation,  respeetively.  The 
amplifieation  eontrol  term,  whieh  modified  amplifieation  rate,  was  modeled  as  a  eombination  of  multiple 
inhibition  terms  and  one  aetivation  term: 


Vamp  —  min  [famp  {TFPI)  ,  famp  (^4)  ,  famp  {^amp))  (18) 

where  Zamp  =  fV  x  fX  x  fVIII  x  fIX.  Although  {Zamp)  is  an  aetivating  term,  we  ineluded 
it  in  the  min  integration  rule;  the  factors  in  Zamp  were  essential  for  amplifieation  (if  any  of  these  factors 
was  missing  the  amplification  reaction  would  not  occur).  Thus,  the  faetors  in  Zamp  were  required 
components,  a  elassification  that  we  implemented  by  the  min  selection  rule.  The  rate  activated  protein  C 
formation  was  given  by: 

Vapc  k  jYp(jfQymation  {TM)  - ^ -  (19) 

^  formation, PC  ‘  X^ 

where  kAPc, formation  and  K formation, pc  denote  the  rate  and  saturation  eonstants  governing  aetivated 
protein  C  formation,  respectively  and  TM  denotes  the  thrombomodulin  abundance.  We  modeled  the 
control  term  which  governed  APC  formation  as  a  single  thrombin-dependent  aetivation  term: 


Vapc  =  max  {f+a  (3^2)) 


(20) 
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Lastly,  we  included  direct  irreversible  inhibition  of  Flla  by  ATTTT: 


(21) 


where  7  was  estimated  to  be  7  =  1.26.  For  ATTTT  inhibition  of  FTTa,  the  control  variables  Vinh,ATiii 
was  taken  to  be  unity.  The  model  equations  were  encoded  using  the  Python  programming  language  and 
solved  using  the  ODEINT  routine  of  the  SciPy  module  [66].  The  model  files  can  be  downloaded  from 
http://www.varnerlab.org. 

4.2.  Estimation  of  Model  Parameters  From  Experimental  Data 

Model  parameters  were  estimated  by  minimizing  the  difference  between  simulations  and 
experimental  thrombin  measurements  (squared  residual): 


(22) 


where  Xj  (r)  denotes  the  measured  value  of  species  j  at  time  r,  Xj  (r,  k)  denotes  the  simulated  value 
for  species  j  at  time  r,  and  coj  (r)  denotes  the  experimental  measurement  variance  for  species  j  at  time 
r.  The  outer  summation  is  with  respect  to  time,  while  the  inner  summation  is  with  respect  to  state. 
We  minimized  the  model  residual  using  Particle  swarm  optimization  (PSO)  [67].  PSO  uses  a  swarming 
metaheuristic  to  explore  parameter  spaces.  A  strength  of  PSO  is  its  ability  to  find  the  global  minimum, 
even  in  the  presence  of  potentially  many  local  minima,  by  communicating  the  local  error  landscape 
experienced  by  each  particle  collectively  to  the  swarm.  Thus,  PSO  acts  both  as  a  local  and  a  global 
search  algorithm.  For  each  iteration,  particles  in  the  swarm  compute  their  local  error  by  evaluating  the 
model  equations  using  their  specific  parameter  vector  realization.  From  each  of  these  local  points,  a 
globally  best  error  is  identified.  Both  the  local  and  global  error  are  then  used  to  update  the  parameter 
estimates  of  each  particle  using  the  rules: 


—  diA.i  +  02ri  (Ci  —  kj)  +  03r2  {Q  —  kj) 
k,  =  kj  +  Ai 


(23) 

(24) 


where  (^i,  62, 6f)  are  adjustable  parameters,  Ci  denotes  the  local  best  solution  found  by  particle  i,  and 


Q  denotes  the  best  solution  found  over  the  entire  population  of  particles.  The  quantities  ri  and  r2 
denote  uniform  random  vectors  with  the  same  dimension  as  the  number  of  unknown  model  parameters 


(/C  X  1).  In  thus  study,  we  used  (6*1,  6*2,  6*3)  =  (1.0, 0.05564,  0.02886).  The  quality  of  parameter  estimates 
was  measured  using  goodness  of  fit  (model  residual).  The  particle  swarm  optimization  routine  was 
implemented  in  the  Python  programming  language.  All  plots  were  made  using  the  Matplotlib  module  of 
Python  [68]. 

4.3.  Global  Sensitivity  Analysis  of  Model  Performance 


We  conducted  a  global  sensitivity  analysis,  using  the  variance-based  method  of  Sobol,  to  estimate 
which  parameters  controlled  the  performance  of  the  reduced  order  model  [69].  We  computed  the  total 
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sensitivity  index  of  eaeh  parameter  relative  to  two  performanee  objeetives,  the  peak  thrombin  time  and 
the  area  under  the  thrombin  eurve  (thrombin  exposure).  We  established  the  sampling  bounds  for  eaeh 
parameter  from  the  minimum  and  maximum  value  of  that  parameter  in  the  parameter  set  ensemble.  We 
used  the  sampling  method  of  Saltelli  et  al.  [70]  to  eompute  a  family  of  N  {2d  +  2)  parameter  sets  whieh 
obeyed  our  parameter  ranges,  where  N  was  the  number  of  trials,  and  d  was  the  number  of  parameters 
in  the  model.  In  our  ease,  N  =  10,000  and  d  =  22,  so  the  total  sensitivity  indiees  were  eomputed  from 
460,000  model  evaluations.  The  varianee-based  sensitivity  analysis  was  eondueted  using  the  SALib 
module  eneoded  in  the  Python  programming  language  [71]. 
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